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Abstract 

Coupled oscillator models where N oscillators are identical and symmetri¬ 
cally coupled to all others with full permutation symmetry are found in a 
variety of applications. Much, but not all, work on phase descriptions of such 
systems consider the special case of pairwise coupling between oscillators. In 
this paper, we show this is restrictive - and we characterise generic multi-way 
interactions between oscillators that are typically present, except at the very 
lowest order near a Hopf bifurcation where the oscillations emerge. We examine 
a network of identical weakly coupled dynamical systems that are close to a 
supercritical Hopf bifurcation by considering two parameters, e (the strength of 
coupling) and A (an unfolding parameter for the Hopf bifurcation). For small 
enough A > 0 there is an attractor that is the product of N stable limit cy¬ 
cles; this persists as a normally hyperbolic invariant torus for sufficiently small 
e > 0. Using equivariant normal form theory, we derive a generic normal form 
for a system of coupled phase oscillators with Sn symmetry. For fixed N and 
taking the limit 0<e<CA<Cl, we show that the attracting dynamics of the 
system on the torus can be well approximated by a coupled phase oscillator 
system that, to lowest order, is the well-known Kuramoto-Sakaguchi system of 
coupled oscillators. The next order of approximation genericlly includes terms 
with up to four interacting phases, regardless of N. Using a normalization that 
maintains nontrivial interactions in the limit A —)■ oo, we show that the addi¬ 
tional terms can lead to new phenomena in terms of coexistence of two-cluster 
states with the same phase difference but different cluster size. 


1 Introduction 

Coupled oscillator models are used in a wide variety of applications. They appear in 
neuroscience for studying neuronal oscillation patterns in the brain (see for example 
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mUEllET]); in chemistry (see for example [251125] ) and in physics (see for example 
PEQ]). A powerful method for understanding the dynamics of coupled oscillators 
comes from noting that N limit cycle oscillators give rise to a normally hyperbolic 
invariant torus that persists for weak enough coupling |3]. In such cases we can 
describe the asymptotic dynamics in terms of just phases. A specihc example of a 
coupled identical phase oscillator system with global (all-to-all) coupling is that of 
Kuramoto [IHI 

d K ^ 

=^ +( 1 - 1 ) 

' k=l 

with hxed natural frequency u, coupling strength A" > 0. Although the original 
work of Kuramoto considered g{ip) = sin((p), a more general “Kuramoto-Sakaguchi” 
coupling (phase interaction) function [23] is 


g{ip) = sin((p + a). 


( 1 . 2 ) 


For the system ( 1.1[1.2 ), the only attractors are full synchrony or full asynchrony, 
depending on the value of the parameter a, while in the special case cos(a) = 0 the 
system is integrable. Many papers have studied the dynamics of this and related 
systems; see for example PED- For this permutation symmetric case of identical 
oscillators, the system above is known to behave in ways that are not generic, even 
accounting for symmetries. There can be a large number of integrals of the motion 
[25] and degenerate bifurcation behaviour on varying a [1]. As pointed out in [T5] . 
for weak linear coupling of nonlinear systems near Hopf bifurcation, one expects to 
have a coupling function g that is smooth and 27r-periodic PP. That is, the generic 
situation is that all Ak will be non-zero in the Fourier expansion 


9(4=) = sin(t:(i3 + xt) 


(1,3) 


k=0 


where the will decay with fc at a rate that will depend on the smoothness of g. 

However, recent work by Rosenblum, Pikovsky and co-workers has highlighted 
that more complex interactions may be present in coupled oscillator systems, and that 
this may lead to new emergent phenomena such as self-organized quasiperiodicity IZH 
on including an additional damped equation, or on including direct dependence of a 
phase shift Xk on an order parameter [H] 122] . 

The current paper considers generic nonlinear, but fully permutation symmetric, 
weak coupling of N identical Hopf bifurcations. We show, by examining a generic 
normal form for equivariant Hopf bifurcation and unfolding parameter A, that the 
system has an attracting invariant torus for 0 < e <C A ^ 1. On this torus, the flow 
can be approximated by (1.1[1.2) at lowest order, but over a longer timescale it can 
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be better approximated by a system of the form 


d e ^ e ^ 

—ipj =n{ip, e) + ^ g2{^k - ^j) + ^ 93{^k + 2v7j 

k=l k,i=l 

N N 

^ Jp 34{‘2(pk - n - ^j) + Jp 95{^k + ^e- ^m- ^j)- 

k,l=l k,i,m=l 

The frequency Q{ip,e) is a symmetric function of the phases that is close to the 
frequency at Hopf bifurcation of the uncoupled system, and we have coupling via 



5'2(</5) =6 cos((p + Xi) + 6 cos(2(p + X 2 ) 

93{.v) =6cos(v3 + X3) 

9i{.^) =^4Cos(v3 + X4) 

9b{9^) =^5 cos((p + Xs) 

where and Xj depend on A. More precisely, they are determined by 

92 {, 9 }) cos((p + X?) + Ml cos((p + xi) + Ml cos(2(p + xi) 
93 M) =Mlcos{ip + xl) 

9iM) =A^] cos(v3 + xi) 

93M) =MlcosM + x\) 


(1.5) 


( 1 . 6 ) 


for some constant coefficients and xl- A more precise statement that includes the 
suppressed higher order terms is given in Theorem 3^ and Corollary |3.3 Most of the 
discussion, apart from Section 5.3, will assume N is hxed, but we assume the given 


normalization in (1.4) of the sums by N, N‘^ or to ensure non-trivial coupling in 


the thermodynamic limit N ^ 00 . 


Including only the very lowest order terms, we will see that (1.4) reduces to (1.1) 


with coupling (1.2): 


d 


e 

N 


N 


X] + x°)- 


(1.7) 


k=l 


with f2, and x? constants. As (1.4) shows, to the next order we only need to consider 
interaction terms of up to four phases. Each of the smooth periodic functions 9kM) 
for k = 1,..., 5 involves only one Fourier mode, except for g 2 which involve two, and 
f2((p,e) is a symmetric function of the phases. Note that (1.4) has normal form 
symmetry 

(cpi ,..., (pTv) Ml + X, ■ ■ ■, + x) 


for any x G T, in addition to the permutation symmetries Sjy- 

The structure of the paper is as follows. In Section we give an outline of the 
normal form theory for S'at equivariant Hopf bifurcation on where Sjv acts nat¬ 
urally by permutation of coordinates. This action decomposes into two irreducible 
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subspaces of complex type, one of dimension one (corresponding to bifurcation to 
in-phase oscillation) and one of dimension N — 1 (corresponding to bifurcation to 
anti-phase oscillation). We include two bifurcation parameters, A determining the 
Hopf bifurcation e representing the strength of coupling, in regimes where both are 

small. 

Section considers a set of coupled systems undergoing a generic supercritical 
Hopf bifurcation in the case of weak coupling 0 < e <C A 1. The main result is 
Theorem 3^ which is proved in Section by performing an explicit reduction of the 
normal form to an invariant Wtorus represented by coupled phase oscillators. Sec¬ 
tion briefly considers a numerical example of the reduction as well as discussions of 
the consequences of the new interaction terms on fully synchronous and two-cluster 
states. The new terms introduce a particular (quadratic) nonlinearity to the equations 
for the phases of two cluster states, and in Theorem 131 we detail a particular new 
phenomenon. Finally, in Section we discuss some implications of this on the dynam¬ 
ics of all-to-all coupled oscillators near Hopf bifurcation, and we relate to other work 
in the literature that considers more general nonlinear coupling between oscillators. 


2 Equivariant Hopf bifurcation with Sn symmetry 


Suppose we have N identical and identically interacting smooth ((7°°) dynamical 
systems on Xj G [d > 2), generated by the following coupled ordinary differential 
equations: 


d 



Hx{xi) + ehx,e{xi;x2, ..., xn) 


d 


dt^^ 


Hx(Xn) 


( 2 . 1 ) 


The “coupling parameter” e G M is such that the system completely decouples for 
e = 0. We also assume that each system undergoes a Hopf bifurcation of x = 0 when 
a “Hopf parameter” A G M passes through zero for e = 0. 

Without loss of generality we assume that the uncoupled system for x G given 
by 

^x = Hx{x), (2.2) 

has a linearly stable hxed point at x = 0 for A < 0 that undergoes Hopf bifurcation 
at A = 0. Without loss of generality, we can assume that DHx{0) has a complex pair 
of eigenvalues 

X±iu 


where ca 7 ^ 0 and all other eigenvalues fi of DHx{0) satisfy Re^jj) < —r < 0. 

We also assume that, without loss of generality, (xi,..., x^v) = 0 is an equilibrium 
for (A,e) in some neighbourhood of ( 0 , 0 ). As we will be interested in attracting 
behaviour near bifurcation we assume that the bifurcation is supercritical, i.e. it 
gives rise to a small amplitude stable limit cycle for A > 0. 
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Note that the Jacobian of (2.1) at (xi,... = 0 will have the form 

^ DHx(0) + edih\^f:{0) e<^2^A,e(0) • • • e<^2^A,e(0) \ 


ed2hx,e{0) 


y ed2hx,e{0) 


-D-ffA(O) + edihx,e{0) 


ed2hx,e{0) 


(2.3) 


ed2hx^e{0) ■ ■ ■ -D-ffA(O) + edihx^ei^) ) 


where dk represents the Jacobian with respect to the kth. argument. 

We assume that the coupling respects the fact that the uncoupled systems can be 
permuted arbitrarily, i.e. that the system is equivariant under the action of Sn on 
^dN permutation 


a{xi, ...,xn) = {x„-i(i), x^-i(N)), (2.4) 

for any {x\, ..., Xn) G and cr G Sn ■ 

Although Hopf bifurcation in the absence of symmetry can generically be reduced 
to a two dimensional centre manifold, this is not the case here - the action of the 
symmetry group S'at means that for e > 0 the centre manifold at generic bifurcation 
will generically be either 2 dimensional or 2N — 2 dimensional. In the uncoupled case 
A = e = 0 the extra structure means that the centre manifold will be 2N dimensional. 


3 Hopf normal form for a weakly coupled system 


Using equivariant bifurcation theory [13] it is possible to write the system of ODEs 
) on a centre manifold (zi,..., zj^) G where in the case A = e = 0 the centre 
manifold in each coordinate x^ is parametrized by Zk- This system on the centre 
manifold is 



= fx{zi) + egx{zi,Z2,...,ZN) + 0{e^) 

dt 


— f\{zN) + ^g\{Z]\i]Zi, . . . ,Zf<!_i) + 0{e^) 


(3.1) 


where z G and we note the right hand sides can be chosen to be with r 
arbitrarily large, in a neighbourhood of the bifurcation. The conditions for Hopf 


bifurcation mean that for (3.1) we have /o(0) = 0 and dfo{0) has a pair of purely 
imaginary eigenvalues izcu that pass transversely through the imaginary axis with 
non-zero speed on changing A. The action of Sn on where a E Sn acts by 
permutation of coordinates 


a{zi, ...,zn) = {z„-i(i),. .., z^-i(N)), (3.2) 

where {zi,, zn) G and so gx{zi] Z 2 ,..., zn) is symmetric under all permutations 
of the last — 1 arguments that £x the hrst. 
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Poincare-Birkhoff normal form theory means that to all polynomial orders we can 
assume there is a normal form symmetry given by the action of on 


9{zi, ^ , 2;y). 


(3,3) 


The symmetries (3.2), (3.3) restrict the possible terms that can appear in the 


normal form; we can characterise these by hnding the possible equivariants, one 
order at a time. This can be expressed in the following form which can be recov¬ 
ered from HD], where denotes denotes and Yi,j,k denotes 

E N 

i=l 2-^j=l 2-^k=l‘ 

Theorem 3.1 Suppose N > 4. Let f : —)■ C'^ be Sn x -equivariant with respect 


to the action (3.2), (3.3) with polynomial components of degree lower or equal than 


3. Then / = (/i, /a,..., /at) where 


11 


fl{zi,Z2, . . . ,^Ar) = ^ aihi{zi,Z2, ...,Zn) 


i=—l 


Mzi, 22, ■ ■ ■ , Zn) =/i( 22, Zl, ■ ■ ■ , Zn) 


(3.4) 


and 


fN{zi, Z2,..., Zn) =/l(^Ar, 2:2, • • • , Zi) 


h-i{z) -- 

N 

ho{ 

= 

Zl, 

hi{z) = 

\zi\^Zi 




h2{z) = 


hsi 

:^) = 


hi{z) = 


h^{ 

= 


he{z) = 


hji 

:^) = 

Zljf Yi,j ZiZj, 

hs{z) = 


hgi 

= 

]V2 Sij ^ki 

hio{z) = 

Yi,k , 

/ill 

(:^) = 

W ZiZjZk 

Oj G C. 

Also we denote p 

= 

'j for j 

= l,...,iV. 


(3.5) 


Proof: For details, see [101 Section 2.1.2]. 


We summarise so far: if system (2.1) is such that (a) the system decouples for 


e = 0 and (b) for e = 0 each system has a generic Hopf bifurcation at A = 0, x = 0, 
then near A = e = x = 0 the dynamics can be written on a centre manifold of 


dimension 2N as (3.1). We now state the main result of our paper: 
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Theorem 3.2 Consider system (3.1) with Sj^-symmetry (for fixed N) such that the 
N uncoupled systems (e = 0) undergo a generic supercritical Hopf bifurcation on X 
passing through 0. There exists Aq > 0 and eo = eo(A) such that for any X G (0 , Aq) 
and |e| < eo(A) the system (3.1) has an attracting -smooth invariant N-dimensional 
torus for arbitrarily large r. 

Moreover, on this invariant torus, the phases (pj of the flow can be expressed as a 
coupled oscillator system 


N 




N 


dt 


k=l 


+ Jp dsiTk + Te - 2(^j) 
k,l=l 


N 




k,e=i 


N 


2 ^ gfi2p>k -Te- Tj) + ^ 9^\Tk + Tj) 


k,l,m=l 


+ ■ ■ ■ 1 Tn) + O(e^) 

for fixed 0 < X < Xq in the limit e —)• 0, where fi((^,e) is independent of j and 

92{t) =^i cos((^ + xi) + A^J cos{p> + xl) + cos( 2 (^ + X2) 

93{t) =A^3 cos(<^ + Xs) 
gfiT) cos{(p + xl) 

9s{t) =Klcos{^ + xD- 

The constants ffi and xl generically non-zero. The error term satisfies 

9(.‘Pu---,‘Pn) 


(3.7) 


uniformly in the phases ipk- The truncation of (3.6) by removing g and O(e^) terms is 
valid over time intervals 0 < t < i where i = 0(e“^A“^) in the limit 0 < e -C A <C 1. 
In particular, for any N, this approximation involves up to four interacting phases. 


The proof of this Theorem is given in the next section. We remark that if we 
set = 0 in the theorem above, this gives the Kuramoto-Sakaguchi system as a 
truncation but with errors O(eA), meaning the timescale of validity of the Kuramoto- 
Sakaguchi system approximation will typically only be 0(e“^A“^). We discuss the 
implications on timescales of validity of the approximation more precisely in the 
following corollary which is obtained by integrating the 0{eX^) error term in the 
truncation. 


Corollary 3.3 Consider the system and hypotheses as in Theorem 3.2 Then for any 
0 < a < 1 and any 0 < e -C A 1 such that there is an attracting N-torus, there 
is a timescale tmax = 0(e““^A“^) such that for any solution (pit) of (3.6), there is a 
solution (pit) of the truncated eguation with 

\(p{t) - (p{t)\ < e“ 


overO < t < tmax- If we truncate further to only Kuramoto-Sakaguchi terms by setting 
= 0 , then this will be possible only over a shorter timescale tmax = 0 (e““^A“^). 
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Proof of Theorem 


3.2 


We write the equation for ^zi from (3.1) in Poincare-Birkhoff normal form [13] as 
the S'tv X S'Bequivariant system 


= U{zi) + eFi{zi,... ,ZN,e), 


(4.1) 


and the equations for the other ^^Zj are obtained by permutation of the indices; there 
is an error term but this is beyond all (polynomial) orders. 


Since we are assuming there is a Hopf bifurcation of (3.1) for e = 0 (the uncoupled 
system) on varying A through 0, it follows that 


d 


= U(zj) := V{zi)zj := [A + iiz + ai\zi\‘‘ + r(ji)] Zi, 


(4.2) 


and we write V{zi) = Vr{zi) + iVi{zi). The tail t{zi)zi represents the higher order 
terms in the normal form for the uncoupled system: we can assume r(0) = r'(0) = 
r"(0) = 0 and write t{zi) = TR{zi)+iTi{zi). The hypothesis that the Hopf bifurcation 
is generic and supercritical implies 


aiR < 0 . 


We seek solutions of (4.2) of the form 


z^{t) = 


(4.3) 


for some i?i(t), t/’i(^) ^md constant H. Substituting this into (4.2), we require 


d 

—+ iRi 

dt 






From (4.2), note that 

Vr{Ri) = a + uirRI + tr{RI), Vi{Ri) = u: + aiiR\ + r/(i?i). 


From this, it is clear that for small enough A > 0 and e = 0 there is a stable periodic 
orbit at fixed = i?* > 0 such that Vr{R^) = 0, with angular frequency H = Vi{R^) 
and arbitrary but fixed phase 'ipi. 

More precisely, solving Vr{R^) = 0, we note (recalling aiR < 0) that 


Ri 


A 


— O^lR 


+ 0(A^), 


an 


(4,4) 


n =V,(Ri) = u + auRi + t{R.) = u + —l^A + 0(A"). 


-0-lR 


In particular there is a Aq > 0 such that for any 0 < A < Aq there is a stable periodic 
orbit (|4.3|) satisfying (|4.4[). 














Now consider the dynamics of the fnll (bnt still nnconpled) system. For e = 0 and 
any 0 < A < Aq there is a stable invariant torns given by 




(4.5) 


parametrized by the phases {'ipi,.. .'i/jn) G This invariant torns is foliated by 
nentrally stable periodic orbits with period 271/Q and so for each 0 < A < Aq, the 
torns is normally hyperbolic. By Fenichel’s theorem [12] there is an eo (depending 
on A) snch that for 0 < e < Cq the invariant torus persists and is C"’-smooth for 
arbitrarily large r. Note that reducing r will restrict the eg: we will need r > 5 for 
the approximation to be valid. 

We now aim to hnd the approximating system ( 3.6[ ) on this invariant torus for 
0 < e < cq. We follow a method similar to |29l Section 7.3], using coordinate changes 
and a slow time to “blow up” the weak hyperbolic dynamics. Including all terms up 
to cubic order (except for the linear term aiZi which can be assumed to be contained 
in U{zi) by a suitable change in parameters), using Theorem 3.1 we have 


= 


1 


'n 






3 

’iV 2 


N 






3 

N 




1 


ATS 


hk 


flio 


iV 2 


E 

j,k 


Zj\Zk\ 


'N 


+ “11^5 


iV 2 


2 :^ 


(4.6) 


E 

3,k/ 


ZjZkZi 


3,k 

+ Fi + 0(e). 


where the e = 0 error term is Fi = 0{\z\^). The complex normal form coefficients 
can be expressed using real quantities and 9k (or qur and aui) such that 

ak = = RkR + ittki, 

for i = —1,1,...,11. We seek solutions of the following form: 

Zk{t) = = [i?, +p,(t)]e'(^*+^'=W) 


that remain close to periodic orbits on the invariant torus (4.5). In particular, we 
seek solutions such that pk is small and V’fc varies slowly with t. Re-writing (4.1), we 
have 


—Pi -|- iRi 

at 




= U(R,)+eFi(zu--- ,2«,0)e-<“+''''>+O(e2). (4.7) 


Writing U in real and imaginary parts and expanding for small pi we have 
U{R,) = UR{R,)+tRMR,) 

= Ur{R^: -|- Pi) -|- iRiVi^R^ + pi) 

= (7)j(R*)pi + iRi[Vi{R,) + V/(R*)pi] + 0{pl). 
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If we define 


(4.8) 


A(A) := B(X) := 


A 


A >/2 ’ 


then, from (4.4), 


U(Ri) = Aa1(A)pi + ifiip + A‘''^B(A)pi| + 0(pl). 


(4.9) 


This implies that (|4.7|) can be expressed as 

= XA{X)pi+iRi[Q +X^/^B{X)pi] 

+eFi(zi,...,z^)e-*(™+^i) + 0(e2) 


—Pi + 

at 




(4.10) 


Recalling from (|g that Rl = A/(-ai^j) + 0(A2), U{R,) = Ur{R,) + iVi{R,)R, = 
(A + aiuRl + t{R^))R^, t{z) = 0{z^), and t'{z) = O(z^) so we have 


71(A) = 


UniR*) 

X 

X + SaijiRl + 't'r{R*)R* + 'Tr{R*) 

A 

l + ^^ + 0 (A) 

— 0>1R 

-2 +0(A) 


Similarly, we have 


B{X) = 


Vji^R^Xj 2R^:aij + Tji^R^Xj 


yx 

2a\i\/~X 

's/XaiR 

2aii 


V~'^iR 

In particular, for A —)■ 0 there are hnite limits 

71(0) =-2, R(0) = 

Note that 

a_ih_i(; 2 )e-*(^‘+^i) = 


yx 

(1 + 0(A)) 

+ 0(A). 


2a 


II 


yf—aiR 


(4.11) 


(4.12) 


(4.13) 


^ ^ .gi(Ot+pi)g-i(Ot+pi) 

= a-iR Rj[cos{9_i + V’i - -01) + sin(6'_i + ipj - ipi)]. 
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Applying similar expansions for the remaining terms in Fi and taking real parts of 
(|4.10[) gives 


(^) = 


AA(A)pi + e a_i Rj cos(0_i + ^ipj - 
+ 0^2 RlRj cos( 6'2 + -01 - i’j) 

+ 0^3 Y'j RlRj cos( 6'3 + ijj - tpi) 

+“4 RlR] cos 6^4 

+«5 Yj,k RlRjRk COs(6'5 + ^jJj - 'Ipk) 

+06 Yj R^R] cos{6q + 2^j - 2V’i) 

+«7 Yi,j RiRiRj cos[6'7 + {'ipi - V’l) + {'ipj - V’l)] 
Attg Yl R] cos( 6'8 + V'i - V’i)+ 

+“9 Yj,k R^Rk cos( 6'9 + 2V’j -i^k- i>i) 

+ttio Y'j,k RjRl cos(6'io + 'Ipj - -01) 

+ “11 Y'i i k RlRjRk COs(6'ii + - 0 * + - 0 fc - 0 l) 

+0(p^e^’)’ 


(4.14) 


where = maXj(p2) and Oj := jj Y^=i Eyfc ^j,k ■= Syfc=i ^^jF etc are the 
normalized sums. The equivalent equation for 0i is obtained by taking imaginary 
parts of (|4.7l): 


+ ^0i(t)] — RiVt + i?iA^'^^i?(A)pi + e \oi-i Yj Rj sin(0_i + 0i — 'ipj) 
+02 Y!j RlRj sin( 6'2 + 01 - tpj) 

+Q !3 Yi RlRj siii( 6'3 + - 0l) 

+“4 RiR^n sin Oa 

+“5 Yj,k RlRjRk sin(6'5 + ijjj - ipk) 

+06 y!. RiR] sin( 6'6 + 2 ( 0 J - 0 i)) 

+“7 Yi,j RiRiRj + (0i - 0i) + {il)j - 0i)] 

+08 Y!j R^j sin( 6'8 + - 0 i) 

+“9 Yj,k R^Rk sin( 6'9 + 20J -i/Jk- 0i) 

+“10 Zip RjRl sin( 6 'io + - 0 i) 

+“11 Yij k RlRjRk sin(6'ii + + ijj - ilk- ill)] 

+0(p^e2) 
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which, after cancellation and dividing throngh by i?i, gives 

= X^/^B{X)pi + e[a_i Y^j{Rj/Ri) sin(6'_i + V’l - i^j) 

+Q!2 RiRj sin(6'2 + '01 - i^j) 

+Q!3 Y.j -Ri-Rj sin(6'3 + - 'i/'i) 

+«4 -R| sin 04 

+“5 fc sin (05 + ijj - ipk) 

+06 Z] ■ Rj sin (06 + ‘ 2 {'ipj - 0 i)) 

+«7 Z^ij RiRj sin [07 + (0i - 0i) + {'ipj - ipi)] 

+«8 HaR^/Ri) sin(08 + - i’l) 

+“9 Y.j,kiR]Rk/Ri) sin(09 + 20J - 0fc - 0i) 

+«io Y!i,ki.RjRl/Ri) sin(0io + - 0i) 

+“11 Y.ijki^iRjRk/Ri) sin(0ii + 0i + -ipk- 0i)] 

+ iO(p^’ 62 ) 

We now define scaled radial variables and a slow time T by 


R*{^) rj. 

Pk = e ^ rfc, T = At. 


For fixed A, note that pk = 0(e) and so (4.14) can be written 


where 

and 


^ri(t) = Al(A)ri + /i + 0(e), 


/i:=/i° + ^^(A)/^ + 0 (A 2 ) 

fl ■= «-l T!j COS(0_1 + 01 - 0j), 

/i^ := a 2 Ej cos (02 + 01 - 0 i) 

+"3 E'- cos(03 + 0J - 01) 

+Q!4 E^' cos 04 

+“5 Eyfc cos(05 + 0J - 0fc) 

+«6 El COS (06 + 2 ( 0 j - 0 i)) 

+“7 Eij COS[07 + (0i - 0i) + (0j - 0i)] 
+08 E^' COS(08 + 0j - 0l) + 

+“9 Eyfc COS(09 + 20J - 0fc - 0i) 

+“10 Ep COS(01O + 0j - 0l) 

+“11 Eij,fc COS(011 + 0i + 0j - 0fc - 0l). 


Similarly, (4.16) can be written 


where 


^^0i(0 — eA ^0(A)ri + eA ^hi + O(e^) 


C(A): = WMF = _2^11 + 0(A), 
V A o,iR 

hi: = h° + i? 2 (A)/,J+ 0 (a 2 ) 


(4.16) 


(4.17) 


(4.18) 


(4.19) 


(4.20) 
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and 


h\ := a_i Y.j sin(6'_i + V'l - 

h\ := 02 Y!j sin( 6'2 + - i^j) 

+03 sin( 6'3 + ‘ipj - - 01 ) 

+«4 sin 64 

+«5 Eifc sin( 6'5 + ^pj - 'ijjk) 

+06 E,- sin( 6'6 + 

+«7 Ei,j sin[6'7 + (^i - ^ 1 ) + ii’j - i>i)] 

+08 Yl'j sin(6'8 + V'i - '00 
+«9 Ej,fc sin(6'9 + 2'ijjj -i)k- 00 
+«io Ep sin(6'io + Oi - 00 

+«ii sin(6'ii + Oi + Oi - Ofc - 00 


In snmmary, we can 


write system ( |4.18[4.20D 


as 


^rj - A{\)rj + fj + 0(e) 

^Oi = ^ + hj] + 0{e^) 


(4.21) 


(4.22) 


for j = 1 ,..., A^. Note that A, C, fj and hj have hnite limits as A —)■ 0 and so (4.22) 
gives a slow timescale for evolntion of Oj as long as 


e = o(A) 


which holds, for example, if e = A^. Defining a new set of amplitnde variables 


aj :=rj + 


/00i,...,0jv_0 

Al(A) 


system (4.22) for fixed A becomes 

= Al(A) 


cTj + 0(e) 


0(A) 

a-^ 

^3 A(A) 

+ hj 


+ O(e 0 . 


which gives 

where 

Let ns write 

where 


— A{X)aj + 0(e) 

= e\-YC{X)a,+Hf+0{Y), 

H -h -^f 
Hj = + XHj + 0(AO, 



-h? - 

3 




A 


Al(0) 


fO 

J 7 ? 


m C{0) , 

^ A ( 0)^0 


0 '( 0 )A( 0 )-AI'( 0 ) 0 ( 0 ) ,0 


Al( 0 )^ 


r 

J 7 ? 


(4.23) 

(4.24) 

(4.25) 

(4.26) 

(4.27) 
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and these are trigonometric polynomials in ~ snch that only involves pair¬ 
wise conpling (and on a_i while Hj inclndes conpling of np to fonr phases (and on 

a2, ..., ttii). 

Applying Fenichel’s theorem [T^ Theorem 9.1] to 
is a 0 < eo = o(A) snch that whenever |e| < Cq the evolntion of the phases is given by 

= eX-^Hj + 0{€^) = €X-^[H^ + XHj] + O(eA). 

The solntions of this rednced eqnation are approximated by solntions of 


(4.24), for all 0 < A < Ao there 


A 


tj 


eX-^H, 


eX-^[H^ + XHj] 


over an interval of time 0 < T < T with T = 0{e ^X ^). In terms of the original 
phases (pj and time t, this rednced eqnation is 


d 


—ifj -Q. + e[Hj XHj] 


dt 


(4.28) 


where the phase differences ipj —'4)k = P>j — for all j and fc, and the approximation 
will be close for times 0 < f < t with t = 0{e~^X~^). 

For k = —1,1,..., 11 we dehne fdk and jk snch that for all 6 

(3k cos( 7 j 6) := ak sin(6'fc 6^) - cos{6k + 6). 


Then we can write (4.28) in the form 
d 


-a + iH, 


N 


^ - p>j) + ^3 A + ^e- 2¥^i) 

k=l k,i=l 


N 


(4.29) 


N 


N 


+ Jp ^ dii^tpk - (pe - (Pj) + ^ g5{^k + n - p>m - ^j) 


k,e=i 


N3 


k,l,m=l 


which is eqnivalent to a trnncation of (3.6) snch that the conpling is as follows: 
i2((p, e) =12 Rle 




(di cos 74 + + V’j - 


j,k 


g2{p>) =( 3-1 cos(7_i + p>) + [/?2 cos(72 - p>) + ( 3 ^ cos(73 p>) 

+(3e cos(76 2p>) + (3s cosijs + l^) + (3io cos(7io + p>)] 

C"( 0 )A( 0 ) - ^'( 0 )^( 0 ) 


A( 0)2 

gsi^’) =Rl AcOs(77 + (p)] 
gi{ip) =Rl [(3g cos(79 -1 - ip)] 

g^^ip) =Rl [(3ii cos(7ii + ip)]. 


-a_i cos(6'_i -|- ip) 


(4.30) 
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Substituting in the leading order of the expansion = A/(—Oi/j) + O(A^) and 
simplifying the trigonometric expressions, we obtain the leading order terms in a 
phase description of the bifurcating solutions, involving the terms as expressed in 


(3.7). Note that the a_i term in g 2 involves (via v4'(0) and C"(0)) coefficients at hfth 


order in the normal form for (4.2) 


5 Examples and consequences 

5.1 A numerical example 


We briefly give a numerical example that illustrates the reduction in Theorem 3.2 
Consider the system of N globally coupled Stuart-Landau oscillators of the form (4.1) 
where we choose all parameters zero except for 


A = 0.1, o; = 1, £ = 0.5, air = —1, 0 ‘ 2 r = 0.3 


(5.1) 


Figure shows some time series for the system (4.1) with these parameters along 
with the case a 2 r = —0.3; (a) shows a stable antiphase solution while (b) shows a 
stable in-phase solution. 

Using the results above we re-write in terms of the phase only equations. In case 
(a) we have i?* = a/OT = 0.3162, a 2 = 0.1, 62 = 0, /32 = 0.2, 72 = —vr/2; all of the gu 
are zero except for g 2 {^) = 0.3 sin(93). In case (b) we similarly get ( 72 ( 7 ^) = —0.3 sin((y9). 
Figure [^presents two time series for the system (1.4)). Observe the qualitatively (and 


indeed quantitatively) similar behaviour of (4.1) and the reduced system (1.4) in this 
case, even though A and e are comparatively large. 


5.2 Non-pairwise coupling and synchrony 

We discuss certain periodic cluster states, and hnd that the presence of the additional 
“non-pairwise coupling” terms in (1.4) implies a wider range of behaviour that is 
possible for pairwise coupling (1.1). Note that particular isotropy subgroups for 
S'tv X symmetry that generally contain periodic solutions |3]. Using the notation 
C = we recall from [5] that in particular the following hxed point subspaces 

are invariant: 

Fix(S'Af) = {{z,... ,z) : z e C}, 

Fix(Z7v) = {{z, zC, zC ,..., zC^) : z eC} 
and indeed for any factorization N = km the following hxed point subspace is also 
invariant: 

FixiiS^r ^Zm) = {{z,...,z, zC\ ..., <-") : e C} 

The above hxed point spaces all contain a periodic orbit after the Hopf bifurcation 
(A > 0), for sufficiently small e, even on inclusion of non-pairwise coupling terms. 


This can be verihed by considering the reduced equations (1.4) and noting that within 


each of these invariant spaces we can reduce to a phase equation $ is a constant that 
depends on the coupling: for example, for full synchrony we have 

Fix(^^) = {(4>,...,4))} 
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and so, using ( 4.29[4.30 ) we can determine that the frequency of the fully synchronous 
state is 


11 


4> = SJ + COS7-, + Ml ft cos( 7 ft - 


i=2 


74(0)^ 


In principle also the stability of the synchronous state can be determined, as can the 
frequencies and stabilities of the other periodic orbits. 


5.3 Non-pairwise coupling and two-cluster states 

We study two cluster states in more detail: these have isotropy subgroup Sq x Sp 
with Q + P = N. In this case we have 

Fix(S'Q X Sp) = {( y^i, ■ ■ -,^ 2 ) ■ ^ 1,^2 e M}. 

Q P 


For simplicity we assume here that t{z) ;= 0 in the uncoupled system (4.2) so that 


74 '( 0 ) = C"(0) = 0. Restricting (4.29) to the two-cluster subspace and suppressing 
the O(e^) terms we get 

01 = ip2, Q, P) 

02 = (p2, Q, P) 

where Hi and H 2 are given in Appendix If we write the phase difference between 
the clusters a.s ipi — ip2 then 


^ = e[Hi{^i, (f2, Q, P) - H2{(pi, (f2, Q, ^^)]. 

This ODE can be written 


T = eG(\I') : = 2 esin — 


^ ^ 3 ^ 3 ^ 

Ai cos-h i 3 i sm-h A 2 cos-h B 2 sm — 

1 2 "^ 2 ^ 2 ^ 2 


(5.2) 


(5.3) 


The coefficients Ai,Bi, i = 1,2 can be expressed in terms of a = (Q — P)/N with 
a G (—1,1) and the (3j, 7 ^ for j = 2, 3, 6 ,..., 11 according to the following: 


Ai = Rl 


3 -|- ^ . 3 — ^ . 1 -|- 


-0iisin7ii 


-02 sin 72 

i= 3 , 6 , 8 , 

a + 3a^ 


3,6,8,10 / 


^ -07 sin 77 + 2 

+ 0 _isin 7 _i 


09 sin 79 


Bi = R 


A 2 = RI 


4 

1 — 


011 cos 7 ii- 1 -a 0 j cos 7 j I -|-q; 0 _i cos 7 _i 

7=2,3,6,7,8,9,10 / 

1 oS 

011 sin 7 ii H-- —07 sin 77 + 06 sin 76 


B2 = Rl 011 cos 7 ii-F Q (06 cos 76 -F 07 cos 77 ) 
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where for details we refer to Appendix These expressions allow us to draw some 
conclusions about two-cluster states and, in particular, the influence of non-pairwise 
interactions. 

• Firstly, we note that the factor of sin(T/2) corresponds to there always be¬ 
ing a solution G(0) = 0: this corresponds the fully synchronous solution with 
symmetry Sn- 

• Secondly, in the special case N even and P = Q {so that a = 0), Bi = B 2 = 0 
meaning that there is also a solution G{7i) = 0: this corresponds to a solution 
with symmetry (S'p)^ ® Z 12 . 

• Thirdly, loss of linear stability of synchrony is associated with change in sign 
of G^O). More precisely, synchrony is stable if Ai -|- A 2 < 0 and unstable if 
Ai -f- A 2 > 0. 


• Fourthly, if there is a root G(\h) with T 7 ^ 0 mod 27r this is a non-trivial two 
cluster state. This bifurcates from the fully synchronous solution where G'(0) = 
0. The only possibility of there being NO non-trivial two-cluster solution is at 
such a bifurcation point, such that G'( 0 ) = 0 and G(T) has the same sign for 
all T G (0,27r). 

Moreover, the dependence on a of A* and B^ means we can conclude the following 
result about the set of possible two-cluster states. Note that the case 


^7 — Pd — Pll — 0 


(5.5) 


corresponds to the case of there being only pairwise coupling in (1.4). 


Theorem 5.1 For system (5.3) in the special case of pairwise coupling (5.5) and any 
To G (0, 27r) either 

• G is independent of a and G(To) = 0 for all a G (—1,1), or 

• G(To) = 0 for at most one a G (—1,1). 


In the more general case where (5.5) does not hold, there can be an additional case 
• To is a root o/G(T) = 0 for two distinct a G (—1,1). 

Moreover, this case does appear for certain choices of parameters. 


Proof: This follows on noting that if there is pairwise coupling then the dependence 
on a is Ai = oi, Bi = abi, A 2 = 02 , B 2 = 062 - This means that non-trivial roots To 
of G(T) = 0, from (5.4) must satisfy 


a 



+ &2 sin 



T 

-Oi cos — + 02 COS 


3T 

~2' 
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and so either Bi = B 2 = Q for all a and d'o is a root of the right hand side, or there is 
precisely one real a that satishes this eqnation. If this qnantity satishes — 1 < a < 1 
then there is a corresponding two clnster state with close to this phase difference, for 
large enongh N (exactly this phase difference if a is rational). 

On the other hand, if ( |5.5 ) is not satished there may be qnadratic dependence of 
G on a leading to the possibility of two a. For example, consider the specihc case 
where the normal form coefficients are snch that Ai = 1/8 + Bi = —3a/4 and 
A 2 = B 2 = 0. This has nontrivial roots T where 


3 . T , 

--asm — + - + a' 

4 2 


T 

cos—. 


namely there are roots at 4/ = 7r/2 when a^ — 3a/4 + l/8 = 0. This implies there are 
two cluster states with this phase difference for the two isolated values 


a e {1/2,1/4}. 


QED 


Theorem |5 .1 1 highlights a particular restriction in the existence of two-cluster states 
that holds for very general systems of the form (1.1). Only on addition of additional 


interaction terms of the form shown in (1.4),(1.5) do we start to hnd the sort of 
behaviour one would expect of a generic symmetric system on T^. 


6 Discussion 


In summary, we study a system of N identical systems near generic Hopf bifurcation 
that are symmetrically and weak coupled. In such a case, reduction to a phase 
description will be possible for some neighbourhoor, but a pairwise-coupling models 


such as (1.1) may miss a number of qualitatively different terms. 

Note ( 1.4[1.5 ) can be approximated by ( 1.1[1.2 ) but this approximation and the 
true solution may move apart over a timescale of order (9(e“^A“^). The next approx¬ 
imation includes two, three and four-phase interactions, from cubic nonlinearities in 
the equations and this will be valid for the longer timescale 0(e“^A“^). These terms 
will be important especially near a secondary bifurcation where well-known degen¬ 
eracies of (1.1[1.2) will be unfolded. 


We cannot guarantee that the truncation (4.29) has the same qualitative dynamics 


as (3.6) unless the dynamics of the former is robust to addition of higher order terms. 


Although A is small, we assume e must be smaller in order for the reduction to hold - in 
the event that this does not hold then there may be solutions where the amplitude of 
the different oscillators may vary considerably, and the picture of possible dynamics 
may be much richer ng. In addition to the restrictions on the number of phases 
interacting, the form of the functions is illuminating - (4.30) reveals that while §2 has 


hrst and second harmonics, g^-g^ have only hrst harmonics to this lowest order. 


The presence of the phase-dependent frequency detuning n((p, e) in (1.4) is some¬ 
what surprising. This term is invariant under permutations of the arguments of (p 
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and so does not give any effect in the phase difference dynamics and hence on the 
synchrony properties of the system. However, it will be a measnrable effect, for ex¬ 
ample, affecting the freqnencies of different states of synchrony that will be missing 


from the system (1.1) 


If we look at the special case of pairwise (or the even more special case of linear) 
conpling of nonlinear systems nndergoing Hopf bifnrcation, of the form 


= f{zi) + e^g{zi,. 


N 


then clearly this will give g^, g^ and g^ identically zero. While the g 2 terms in (1.4) 
may be qnite complex and can have higher Fonrier modes - previous work on this 
pairwise coupling has demonstrated present of a rich range of behaviours including 
robust attracting heteroclinic cycles (slow switching) [3] and arbitrary cluster states 
pun]. We remark that terms present in g 2 are also present in the work of Hansel et 
al [12] and the generalization [2|. 


The example we give in Theorem |5.1| is a new but subtle dynamical effect that can 
appear for non-pairwise coupling. We expect there are more remarkable implications 
of non-pairwise coupling waiting to be discovered. For example, there are systems of 


pairwise coupled system of the form (1.1) with = 4 or more identical oscillators 


that possess chaotic attracting states, at least for g{ip) with least four harmonics 
[6]. It is still unknown whether there is a g 2 that gives chaotic attractors for all 
sufficiently large N. We speculate that the additional terms g^-g^ may give broader 
regions of existence of chaotic attractors for ( 1.4[1.5 ) using coupling functions with 
fewer harmonics. 

We hnish by mentioning a couple of examples from the literature that have non- 
pairwise coupling. In the authors consider coupled phase oscillators of the 

form 


N 


dt 


where the order parameter is 


d K ^—V , . 

= a; - — sin((y9i - + aR) 


( 6 . 1 ) 


k=l 


R = 


N 


N 

E 

i=i 


exp rpi 


Writing (6.1) in terms of phases only we hnd 


d 


K 


N 


—(pi = w - — ^ sin I (pi - (pfc -F a 


N 


k=l 


\ 


N 


jp E - v>id) 


( 6 . 2 ) 




which includes non-pairwise terms that include all N phases. This can be derived 
by assuming that the coupling is via a dynamic mean field - note that this justihcation 
assumes there is an extra active degree of freedom in the coupling. In this paper we 
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assume the only dynamic variables are those undergoing Hopf bifurcation. Close in 
spirit to our paper is the analysis of im who consider three communities of oscillators 
with frequencies a;*, i = 1,2,3 such that oji + UJ 2 ~ oj^- They hnd three-phase 
interactions in the phases of the order parameters, using an Ott-Antonsen reduction 
|2Uj ■ By comparison, we are dealing simply with the oscillator phases, and the general 
form of our equations means that we are unable to apply the method of |20j . 
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A Coefficients for two cluster states 


Note that if we write p = P/N and q 
have 


Q/N then (4.30) implies that for (5.2) we 


Hi{q:>i,q:>2,Q,P) =/ 3 -i[gcos7_i + pcos(7_i - (pi + ^2)] 

+ Rl [/?2[gcos(72) +pcos(72 + <Pl - (P2)] 

+ /33[gcos(73) +pcos(73 + (p2 - q^i)] + (^^{p + q) cos 74 
+ + g^) cos75 +pgcos(75 + ((p2 - Pi)) +pgcos(75 + (pi - P2))] 

+ [geos 76 +pcos (76 + 2(932 - <Pi))] 

+ ( 3 -r[q^ cos(77) + 2 pq cos(77 + (932 - Pi)) + P^ cos(77 + 2(932 - 9^1))] 

+ / 38 (gcos 78 +pcos (78 + P 2 - 9^1)) 

+ /39[g^ cos79 +pg cos(79 + 2(932 - 951)) 

+ pgcos (79 - 932 + 931) + p ^ cos (79 + 932 - 934)] 

+ /3io(gcos7io +pcos(7io + 9^2 - ^^i)) 

+ / 3 ii[g^ COS711 + 2 pq^ cos(7ii + 9^2 - 9^1) + QP^ cos( 7 ii + 2(932 - 934)) 
+pg^ cos(744 + Pi - P2) + 2 qp^ cos 744 + p^ cos(744 + 9^2-954)] • 

and H2{pi^ 9^2, Q, P) = Hi{p2, Pi, P, Q) and so if we set 4 ' := 934 — 932 then 

= e[H,{p^,p 2 , Q, P) - H 2 {pi, P 2 , Q, P)] =: eG{^). (A.l) 

This can be written 

G(\l/) =(]_i[{q - p) cos7-4 +pcos(7_4 - T) - gcos(7_4 + T)] 

+ Rl [/92[(g -p) cos72 +pcos(72 + T) - gcos(72 - T)] 

+ /^3[(g -p) COS73 +pcos(73 - T) - gcos(73 + T)] 

+ PgHq-p) cos 76 +pcos(76 - 2 T) - gcos(76 + 2 T)] 

+ ~ P^) cos 77 + 2pgcos(77 — T) — 2pgcos(77 + T)] 

+ / 37 [p^ cos(77 — 24 /) — g^ cos(77 + 24 /)] 

+ (^s[{q-p) COS78 +pcos(78 - T) - gcos(78 + T)] (A. 2 ) 

+ /59[(g^ -p^)cos79 +p^cos( 79 - T) - g2 cos(79 + T)] 

+ Ao[(g -p) cos 740 +Pcos(740 - 4 ')) - gcos(74o + 4 '))] 

+ / 5 ii[(g^ - 2pg^ + 2 qp^ - p^) cos 744] 

+ I3ii[{2pq^ - gp^ +p^) cos(744 - T)] 

+ Ai[(-g^ +Pg^ - 2 gp^)cos (744 + d')] 

+/?ii[gp^cos (744 - 2T) - pg^ cos (744 + 2T)]] . 
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Now writing 1 - cos \1/ = 2 sin^(\I//2), sin\l> = 2sin(\I//2) cos(\l//2) and 1 - cos2\l> = 
2 sin^ we obtain expressions 


: = Aisintl/ + Biil — cosd') 

+ 742(sin2\I/ — sin 'b) + i? 2 (cos 'b — cos2\I/) 

d' r ^ 'b d' ^ dd' 3 '^ 

= 2 sin — Ai cos-h -Bi sin- \- A 2 cos-h B 2 sin — 


(A.3) 


where Ai, Bi, i = 1,2 depend on p,q and the jSj, 7^ for j = —1, 2, 3, 6, ..., 11 as 
follows: 

Ai := /?_i(p + g) sin7_i + / 2 ^[-/ 32 (p + g)sin72 

+ + ‘^P^Q + + q^) sin 711 + /97(p^ + 4 pg + q^) sin 77 

+ I^ 9 {P^ + (f) sin 79 + /3io(p + q) sin710 + /32(p + q) sin 72 
+ I 3 ^{P + q) sin73 + I 3 q{p + q) sin76 + I 3 ^{p + q) sin73], 

Bi := ( 3 -i{q - p) cos7_i + i?^[/?2(g - p) cos72 

+ + ‘^p^q - ‘^pq^ + q^) cos711 + | 3 ^(pf' -p^) cos77 

+ ( 3 g{q^ - p^) cos79 + / 3 io{q - p) cos710 + ( 32 {q- p) cos72 

+ /^ 3 (?-p) COS 73 +/?6(g-p) COS 76 +/ 38 (g-p) COS 78 ], 

A2 ■■= B^l^nip^q + pq^) sin7ii + ( 3 ^{p^ + q^) sin77 + /?6(p + q) sin75], 

B2 := Rl[l 3 ii{pq^ - p^q) cos711 + | 3 ^{q^ - p^) cos77 + jiaiq - p) cos76]. 


Finally, we dehne a = q — p with a G (—1,1) so that p = (1 + a)/2, g = (1 — a)/2 
which gives the expressions in (5.4). 
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Figure 1: Example time series for coupled oscillators (4.1) showing the real parts 
Zkrit) against t for TV = 3 and two different parameter sets; (a) shows evolution to 
an anti-phase solution for parameters (5.1) while (b) shows evolution to an in-phase 
solution for the same parameters except a 2 r = —0.3. 




Figure 2: Example time series for a phase approximation (1.4) of the system in 
FigureIt shows i?* cos((pfc(t)) against t for = 3 and two different parameter sets; 
observed that (a) and (b) show similar qualitative and quantitative dynamics to the 
corresponding plots in Figure 
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